% code to compute MT trfn from raw ts files !!!
% 
% clears unneceesary vriables after
% invoking readrts reads the header
% compute cal. table according to gain settings etc;
% 11.7.01 computes cross and auto spectra & average it 
% 12.7 - MT app rest & phse are now computed

% - bug - phase distortion for 2nd and 3rd band- 05.07.2002 corected 30.10.2002 (hurray !)
% - OK for all bands chacked 1.11.2002
% latest date 1.11.2002
clear;
readrts;

Notchi = NotchF(2:5);		% NOTE ! there may be problems 
IsNotch = Filt(2:5);		%        if notch do not have right dimension..
hpb1 = Filt(1);			% NOTE ! high pass filter is always disabled
prob = [Exd Eyd Eyd Eyd Eyd];	% 	 in calib05.m
gain = [Exa Eya Hxa Hya Hza];
clear Bden Record3 Bdend Bddm Count Record4 hpf Bdst Record5;
clear Bdstart MaxRecord ans nt Filt NBlocks array nlines Events ;
clear NotchF dummy Record1 factor Channel;
MaxVal      = 1.0E+6;
MinVal      = 1.0E-6;
ScaleVal    = 1.0E+4;



block = 1024;

notch = Notchi(find(IsNotch>0)); 
PAI = 3.1415926535897;

for i = 1:nchannel,
[cplx,frq1] = calib05(i,nband,block,prob(i),notch,hpb1,gain(i));
CalData(i,:) = cplx;
end;

%------------native declarations---------------

WI = hanning(block);
len = length(matrix)/block;
start = 0;
%len = 17; % delete it !!!! 2/9/2
st = start*block + 1;
en = (start+1)*block;
nchannel = 5;

[TLine,TLRad,TLFreq,TBand] = mtfreq(block); % target freq, line & radius for all bands

TLine1  = TLine(find(TBand==nband));
TLRad1  = TLRad(find(TBand==nband));
TLFreq1 = TLFreq(find(TBand==nband));

nfrq = length(TLine1);
HxHx = zeros([nfrq,1]);

HyHy = HxHx;
HxHy = HxHx;
HyHx = HxHx;
ExHx = HxHx;
EyHx = HxHx;
ExHy = HxHx;
EyHy = HxHx;


for i = (start+1):len,

% FFT, calibration etc.

for u = 1:5,
	K = CorTre(matrix(st:en,u),0);
	Channel1(:,u) = fft(K.*WI);
	Channel(:,u) = Channel1(1:block/2,u).*conj(CalData(u,1:block/2)');
end; 	



% auto and cross spectra---------------------------
for ii = 1:nfrq,
par = parzen(TLRad1(ii));
HxHx(ii) = HxHx(ii) + doparzen(Channel(:,3),conj(Channel(:,3)),TLine1(ii),TLRad1(ii),par);
HyHy(ii) = HyHy(ii) + doparzen(Channel(:,4),conj(Channel(:,4)),TLine1(ii),TLRad1(ii),par);
HxHy(ii) = HxHy(ii) + doparzen(Channel(:,3),conj(Channel(:,4)),TLine1(ii),TLRad1(ii),par);
HyHx(ii) = HyHx(ii) + doparzen(Channel(:,4),conj(Channel(:,3)),TLine1(ii),TLRad1(ii),par);
ExHx(ii) = ExHx(ii) + doparzen(Channel(:,1),conj(Channel(:,3)),TLine1(ii),TLRad1(ii),par);
EyHx(ii) = EyHx(ii) + doparzen(Channel(:,2),conj(Channel(:,3)),TLine1(ii),TLRad1(ii),par);
ExHy(ii) = ExHy(ii) + doparzen(Channel(:,1),conj(Channel(:,4)),TLine1(ii),TLRad1(ii),par);
EyHy(ii) = EyHy(ii) + doparzen(Channel(:,2),conj(Channel(:,4)),TLine1(ii),TLRad1(ii),par);
end;

%------------------------------------------------

%clear Channel;

st = en + 1;
en = (i+1)*block;
end;



Zxx = HyHy.*ExHx - HyHx.*ExHy;
Zyy = HxHx.*EyHy - HxHy.*EyHx;

Zxy = HxHx.*ExHy - HxHy.*ExHx;
Zyx = HyHy.*EyHx - HyHx.*EyHy;
Den = HxHx.*HyHy - HxHy.*HyHx;



frq = frq1;


for i = 1:nfrq,
   if abs(Den(i)>0)&abs(Zxy(i))>0&abs(Zyx(i))>0,
      ZXY(i)=Zxy(i)/Den(i);
      ZYX(i)=Zyx(i)/Den(i);
      ZXX(i)=Zxx(i)/Den(i);
      ZYY(i)=Zyy(i)/Den(i);
   end;
end;




for i = 1:nfrq,
   Phxy(i) = atan(imag(ZXY(i))/real(ZXY(i))) * 180 / 3.1415926535897;
   Rhxy(i) = 0.2 * (1/TLFreq1(i)) * (ZXY(i) * conj(ZXY(i)));
   Phyx(i) = atan(imag(ZYX(i))/real(ZYX(i))) * 180 / 3.1415926535897;
   Rhyx(i) = 0.2 * (1/TLFreq1(i)) * (ZYX(i) * conj(ZYX(i)));
end;


%--------------------------------------------------------------
T = 1./TLFreq1;

subplot(2,1,1);
loglog(T,Rhxy,'b.',T,Rhyx,'r.');
ylabel('ohm-m');
xlabel('period (s)');
grid;
subplot(2,1,2);
semilogx(T,Phxy,'b.',T,Phyx,'r.');
xlabel('period (s)');
ylabel('degrees');
title(file);
grid;




